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The thermal conductance in graphene nanoribbon with a vacancy or silicon point defect (sub- 
stitution of C by Si atom) is investigated by non-equilibrium Green's function (NEGF) formalism 
combined with first-principle calculations density-functional theory with local density approxima- 
tion. An efficient correction to the force constant matrix is presented to solve the confiict between 
the long-range character of the ab initio approach and the first-nearest-neighboring character of the 
NEGF scheme. In nanoribbon with a vacancy defect, the thermal conductance is very sensitive 
to the position of the vacancy defect. A vacancy defect situated at the center of the nanoribbon 
generates a saddle-like surface, which greatly reduces the thermal conductance by strong scattering 
to all phonon modes; while an edge vacancy defect only results in a further reconstruction of the 
edge and slightly reduces the thermal conductance. For the Si defect, the position of the defect 
plays no role for the value of the thermal conductance, since the defective region is limited within 
a narrow area around the defect center. 

PACS numbers; 68.65.-k, 61.46.-w, 65.80.-g, 71.15.Mb 
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I. INTRODUCTION 



Defect engineering is a powerful approach to generate 
nano-materials with plenty of novel properties''^. For in- 
stance, the defect induced magnetism in graphene was 
predicted^''^, and a patterned defect can be used to open 
up a band gap in graphene by breaking the sublat- 
tice symmetry of the honeycomb structure^. The de- 
fects formation in graphene can be observed directly 
by high-resolution transmission electron microscopy 
(TEM)^. The formation and anneahng of Stone- Wales 
and other defects were observed by transmission elec- 
tron aberration-corrected microscope monochromated 
(TEAM 0.5) within 1 A resolution at an acceleration 
vohage of only 80 kV that is low enough to protect 
graphene samples from being destroyed^. These tech- 
niques for observation of defects in atomic-scale serve as 
a fundamental support for the nanoengineering of de- 
fects in graphene. The engineering of defects on the 
two-dimensional graphene sheet has been proposed the- 
oretically in 2008 by Lusk et al. using first-principle 
calculations^. New carbon allotropes can be fabricated 
by periodically arrangement of basic defect building 
blocks such as Stone- Wales defects, inverse Stone- Wales 
defects, vacancy defects, and eto^. The inverse Stone- 
Thrower- Wales defect was also suggested as a build- 
ing block for monolithic nanoengineering on graphene^. 
Quite recently, a one-dimensional topological defect in 
graphene is realized in the experiment^. This extended 
defect exhibits similar property as a metallic wire and 
may form building blocks for various one-dimensional 
nano devices. In present paper, we will focus on the 
vacancy defect and substitutional Si defect in graphene 
nanoribbon. 



In the quasi-one-dimensional nanoribbon structure, 
the two edges undergo a reconstruction and play an im- 
portant role in determining many properties'*^' The 
TEAM 0.5 design was implemented to observe the edge 
reconstruction and capture the complex behavior of a 
atom at the edge of the graphene nanoribbon'^. As the 
edge region of a nanoribbon has very different proper- 
ties from the inner region, an edge defect should behav- 
ior quite differently from a defect inner the nanoribbon. 
Particularly, we are interested in how the edge and in- 
ner point defects behavior distinctly on the phonon ther- 
mal conductance in the graphene nanoribbon. The ther- 
mal transport capability is important in many functional 
nano devices. High thermal conductivity is a key prop- 
erty for high performance electronic devices, while ex- 
tremely low thermal conductivity is crucial for thermo- 
electric materials. For graphene sheet, a high thermal 
conductivity was discovered experimentally'^. Various 
theoretical approaches have been applied to study the 
thermal transport in graphene system. The direct non- 
equilibrium molecular dynamics simulation can be used 
to mimic the experimental conditions to study the ther- 
mal transport in graphenei^^— . The Bolzmann equa- 
tion method is suitable for the study of the diffusive 
thermal transport. The non-equilibrium Green's func- 
tion (NEGF) approachi^Ti^S and a molecular dynamics 
simulation with quantum heat baths^' are both proper 
methods for the study of quantum thermal conductance. 
In our calculation, we apply the NEGF approach. 

In this paper, we study the phonon thermal conduc- 
tance in graphene nanoribbon by using the NEGF ap- 
proach combined with the first-principle calculation. The 
ab initio calculation gives a force constant matrix with 
long-range interaction while the NEGF approach requires 
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first-nearest-neighboring (FNN) interaction in the calcu- 
lation of the surface Green's function. This conflict may 
be improved by using an enlarged unit cell. However, the 
enlargement of unit cell distinctly increases the simula- 
tion cost, which is even worse in ab initio calculation. We 
thus propose an efRcient correction to the force constant 
matrix to effectively renormalize the long-range interac- 
tion components into the FNN part. We then investigate 
vacancy and Si point defects in the nanoribbon. The po- 
sition of the vacancy defect plays an important role. An 
edge vacancy defect results in a further reconstruction of 
the edge of the nanoribbon thus only slightly affects the 
thermal conductance. The inner vacancy defect causes 
the defective area to become saddle-shaped, which con- 
siderably reduces the thermal conductance. All phonon 
modes are strongly scattered by this saddle-like surface 
in the nanoribbon. In the case of Si defect, the thermal 
conductance shows no difference for edge or inner point 
defects. We also study phonon modes in the nanorib- 
bon and find that there are some phonon modes local- 
ized at the edge of the nanoribbon. The number of edge 
modes is determined by the number of C-C bonds at the 
edge of the nanoribbon. These edge modes are of much 
higher frequency, since the edge C-C bonds are shorter 
and stronger after reconstruction during structure relax- 
ation. The edge vacancy or Si defects cause a significant 
downshift of a particular edge mode localized around this 
defective area. 



II. COMBINATION OF SIESTA AND NEGF 

In this section, we concentrate on the combination of 
the NEGF approach and the first-principle calculation. 
This is nontrivial, because the first-principle calculation 
gives a force constant matrix with long-range interaction, 
while what the NEGF requires is a force constant matrix 
of FNN interaction. We have to take care of the conflict 
between the long-range and FNN characters. 

In the ab initio calculation, we use the SIESTA 
package^^ to optimize the structure of the graphene 
nanoribbon and extract the force constant matrix. The 
local density approximation is applied to account for 
the exchange-correlation function with Ceperley- Alder 
parametrization23, and the double-^ basis set orbital (DZ) 
is adopted. During the conjugate-gradients optimization 
of the nanoribbon, the maximum force on each atom is 
smaller than 0.01 eV/A. A mesh cut off of 320 Ry is used. 
Periodic boundary condition is applied in the heat trans- 
port direction, while the nanoribbon is free in the two 
vertical directions by introducing large vacancy space. 
We employ Gamma point k sampling in the periodic di- 
rection. 

As discussed in previous works^i, the original force 
constant matrix Kq from SIESTA violates the rigid trans- 
lational symmetry (RTS) and rigid rotational symme- 
try for the lattice structure. RTS is the origin of the 
zero frequency for three acoustic phonon modes at F 



point, while the rigid rotational symmetry is the root for 
the zero frequency rotational phonon modes in cylinder- 
like structures such as carbon nanotubes. In the pla- 
nar graphene sheet, it is important to ensure that the 
force constant matrix satisfles the RTS, so that the three 
acoustic phonon modes have zero frequency. The RTS 
can be enforced by a renormalization of ifcr^- We refer 
to this new force constant matrix as K. The K matrix 
can be written into following explicit form: 
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where we have used labels L/C/R to denote three differ- 
ent regions: left lead, center, and right lead. The division 
of a whole system into three different regions is conven- 
tional in NEGF approacbi^r^. The two matrices at the 
right top and left bottom corners are resulting from the 
periodic boundary condition in the heat transport direc- 
tion. V^*^ and V^'^' are coupling between leads and cen- 
ter. The force constant matrix for left/right leads have 
been written in a periodic form, since the leads are the 
repetetion of the smallest unit cell. The dimension of 
Koo, Koi and Kio is the total degree of freedom in the 
unit cell. Koo is the force constant matrix for each unit 
cell and Koi is the coupling between two neighboring unit 
cells. Essentially, they are used to reconstruct the force 
constant matrix of the leads in the NEGF scheme: 
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(2) 



We show the result of left lead and the right lead is 
analogous. The most important physical quantity in the 
NEGF approach is the surface Green's function of the 
lead, which is determined by the K^. The K^ is & semi- 
infinite matrix as the lead is semi- infinite. As a result, 
it is impossible to calculate the surface Green's function 
from an inverse of the force constant matirx by the ex- 
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pression ql/r 



An efficient it- 



[uj + ivi) -Kl/h 

erative method has been developed by Sancho et al to 
solve the surface Green's function numerically^^. This 
approach takes the advantage of the periodicity of the 
lead. It also requires that each unit cell only interacts 
with its two FNN unit cells. We refer to this requirement 
as the FNN property of the NEGF approach. After ob- 
taining the gL/Ri the retarded Green's function for the 
center region can be calculated by: 



C = [{uj + i7^Y -Kc-^l-^r 
where the self-energy for the lead is 



(3) 



(4) 
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FIG. 1: (Color online) Optimized configurations for graphene 
nanoribbon. (a). Pristine nanoribbon. The unit cell is high- 
lighted by a thick box (red online), (b). Vacancy defect at 
edge, (c), (d) Vacancy defects inner nanoribbon. (e). Si de- 
fect at the edge. (f). Si defect inner nanoribbon. Two arrows 
in (a) denotes the position for vacancy or Si defects. The 
defective regions are highlighted by a thin box (green online). 



The transmission for phonon is given by the Caroli for- 
mula: 



T[oj] = TriG^rLC^TR), 



(5) 



where G" = (C)^ is the advanced Green's function and 
Ta = — Finally, the phonon thermal conduc- 
tance can be obtained by Landau formula: 



^ J dujUujTpklu 



dn{u},T) 



dT 



(6) 



Now we keep track of the RTS in the whole NEGF 
scheme as shown above. Overall, there are two direct 
consequences from the RTS. Firstly, the diagonalization 
of K should give three zero-frequency acoustic phonon 
modes. Secondly, corresponding to these three modes, 
the transmission function from NEGF should be 3 at w = 
0, i.e T[uj = 0] = 3. In the above, as the K matrix in 
Eq. ([T]) satisfies the RTS, we can correctly achieve zero 
frequency for the three acoustic phonon modes. So the 
first consequence of the RTS is obtained. We now check 
the second consequence. As can be seen from the first 
line of K, the RTS results in following relation 



{Koo + Koi+K, + KiQ)Idof = 0, 



(7) 



where Idof is a unit row vector with dimension as the 
total degree of freedom in unit cell. Ks corresponds to 
the long-range interaction beyond the FNN. The nonzero 
Kg is important to keep the whole K matrix satisfying 
RTS. Because of Ks ^ 0, it can be seen from the first line 
of Kl in Eq. ([2]) that RTS has been broken in the lead. 
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FIG. 2: (Color online) The transmission function with and 
without correction to the force constant matrix for pure 
graphene nanoribbon. 



as {KoQ + Koi + Kiq) Idof 7^ 0. It leads to T[uj = 0] ^ 3 
as shown in Fig. [Hby the solid line (blue online). This 
is the conflict between the long-range character of the 
first-principle calculation and the FNN property of the 
NEGF. 

A straightforward technique to improve this confiict 
is to enlarge the unit cell of the leads to obtain smaller 
value of Kg^^. One can get Kg = for large enough 
unit cell. However, larger unit cell leads to the increase 
of the size of whole system, thus the calculation cost is 
increased and even more difficult in the first-principle 
calculation^^. Another possible way to solve the conflict 
is doing correction to K matrix, while the unit cell kept 
unchanged. There is no increase of calculation cost in this 
method since the size of the unit cell is the same. We can 
use the smallest unit cell in the nanoribbon as shown in 
Fig. [T] (a). This correction is realized in following two 
steps. 

(a) . As discussed above, we enforce the RTS to the 
force constant matrix. This will give us the K matrix 
with translational symmetry fulfllled. We stress that K 
has a nonzero matrix Kg. 

(b) . We set Ks = in matrix, resulting in a new 
force constant matrix K' . 

Obviously, K' from step (b) violates the RTS. We need 
to go back to step (a) to enforce RTS to K' matrix and 
will obtain a new force constant matrix K, which satisfies 
the RTS and has a new nonzero Kg ■ This nonzero matrix 
Kg is smaller than the original one. After many cycles 
of these two steps, we obtain a force constant matrix K 
with RTS fulfilled and Kg w 0. So the FNN property 
of the NEGF can also be satisfied. Essentially, the long- 
range component. Kg, is renormalized into the FNN part 
by the correction. As shown in Fig. [21 by dot line (green 
online), after 20 cycles of steps (a) and (b), the trans- 
mission function from NEGF gives a correct result, i.e 
T[ll! = 0] w 3. The correction has a nice property that it 
only affects the low-frequency region of the transmission 
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FIG. 3: (Color online) The convergence of the correction to 
force constant matrix, (a). Frequency of three acoustic modes 
v.s correction steps Nc- (b). The transmisstion function at 
CO = from NEGF approach v.s correction steps. 

function, indicating that the RTS is enforced without af- 
fecting higher frequency region. The convergence of this 
correction is demonstrated in Fig. [3] It shows that af- 
ter 10 cycles, zero frequency is achieved for the three 
acoustic phonon modes from K and the T[llj = 0] « 3 
from NEGF approach. The correction to the K matrix 
can successfully combine the long-range SIESTA with the 
FNN NEGF sheme. 

Before continuation of further numerical calculation, 
we give a brief summation for this part. (1). We are 
discussing the conflict between the long-range character 
of SIESTA and the FNN property of NEGF. The con- 
flict is explicitly manifested by the RTS in the graphene 
nanoribbon. (2). We solve this conflict by effectively 
renormalizing the long-range component, Ks, into other 
FNN components in the K matrix. 



III. CALCULATION RESULTS AND 
DISCUSSION 

We first study the structure relaxation of the defect- 
free graphene nanoribbon as shown in Fig. [T] (a). We 



FIG. 4: (Color online) (a). The transmission function in 
graphene nanoribbon with different defects, (b) . The phonon 
thermal conductance v.s temperature. 



apply periodic boundary condition in the horizontal di- 
rection and the other two directions are free by intro- 
ducing a vacancy space of 15 A. A larger vacancy space 
of 50 A has also been tried, and no obvious difference 
was observed. Reconstruction happens at two edges of 
the nanoribbon. Five different C-C bond lengths are 
obtained in the nanoribbon. Those ten C-C bonds at 
the two edges are about 1.23 A which is about 20% 
shorter than the other inner bonds. So the edge bonds 
are much stronger. For the three bonds around inner 
carbon atoms, the horizontal one is about 1.41 A while 
the other two bonds are about 1.48 A. The transmis- 
sion function for the defect-free graphene nanoribbon is 
shown in Fig. 2] (a) by solid line (blue online) . After do- 
ing correction to the K matrix as discussed above, we can 
obtain correct value oi T[uj = 0] = 3 in low frequency re- 
gion. There are some phonon modes with extremely high 
frequency (around 1900 cm^^), which does not exist in 
two-dimensional graphene sheet. They are characteristic 
for the graphene nanoribbon with two edges compared 
with the two-dimensional graphene sheet. These modes 
are localized at the edges of the nanoribbon. Fig. [5] 
displays the vibrational property of two selective edge 
modes. Edge modes are originating from the stretching 
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FIG. 5: (Color online) The vibrational morphology of two 
edge modes in pure graphene nanoribbon. 



of edge bonds. The extremely high frequency of the edge 
mode is due to much stronger edge bonds compared with 
other inner bonds. There are totally ten edge modes 
in this particular graphene nanoribbon, corresponding to 
those ten edge C-C bonds. The thermal conductance 
from NEGF is shown in Fig. |4] (b), where the value of 
thermal conductance at room temperature is about 1.5 
nW/K. It yields a room temperature thermal conductiv- 
ity of about 2800 W/m/K, taking 3.35 A as the thick- 
ness of graphene nanoribbon and 775 nm as the phonon 
mean free patbSI. This value is close to the first-principle 
calculation reuslt^* and that from a realistic Gruneisen 
parameter—. 

We now study the graphene nanoribbon with one 
vacancy defect situated at two nonequivalent positions 
(edge or center) as indicated by arrows in Fig. [1] (a). 
These two highlighted carbon atoms can be removed to 
generate a vacancy defect in the nanoribbon. Defective 
regions in all structures in Fig. [1] have been highlighted 
by a box (green online). Fig. [T] (b) shows the optimized 
structure for nanoribbon with an edge vacancy defect. 
Reconstruction at the edge is realized by forming a pen- 
tagon, and only a limited region is affected. As a result, 
the edge vacancy defect can only introduce limited effect 
on the transmission function and thermal conductance 
shown in Fig. U) The reconstruction of the edge of the 
ribbon has direct effect on the phonon edge modes. It 
causes a down-shift of about 200 cm~^ for the edge mode 
localized at this defective region. Fig. [T] (c) and (d) ex- 
hibit the relaxed configuration for the nanoribbon with 
a vacancy defect inner the nanoribbon. Two pentagon- 



hexagon 5-6 pairs have been formed around the defect. 
The inner vacancy defect has no effect on the phonon 
edge modes, as the reconstruction only takes place in- 
ner the nanoribbon. Interestingly enough, a saddle-like 
surface is generated around the defective region with the 
defective center as the saddle point. A curvature surface 
was quite common in defective graphene^. The crystal- 
orbital overlap population (COOP)'^" between the defec- 
tive central carbon atom and its four neighboring atoms 
is shown in Fig.|n](b). The COOP for the bonding around 
the vacancy defect is sharper and narrower than the SP2 
bonding in a defect-free nanoribbon as shown in panel 
(a) in the same figure. The bond length of the four de- 
fective bonds is about 1.64 A, which is even larger than 
the C-C bond length of 1.54 A in SP3 bonding diamond. 
The saddle-like surface can efficiently scatter all phonon 
modes in whole frequency region as shown in Fig.|4](a) by 
dash dot line (black online). The high-frequency phonon 
modes are scattered more seriously by this saddle surface 
compared with other types of defects. Only in this inner 
vacancy defect, the low frequency phonon mode will be 
considerably scattered. As a result of the strong scat- 
tering of all phonon modes by this saddle-like surface, 
the phonon thermal conductance is greatly reduced by 
the inner vacancy defect as shown in Fig. 0] (b) by dot 
dash line (black online) . The thermal conductance can be 
reduced by as much as 60% at high temperatures. Con- 
siderable reduction is also observed in low-temperature 
region, since low-frequency phonon modes are also scat- 
tered by the inner vacancy defect. 

After investigation of the vacancy defect, we now 
switch to the discussion of substitutional silicon defect 
in nanoribbon. We consider the two different positions 
(edge/center) as indicated by the arrows in Fig. [1] (a). 
Fig. [T] (e) shows the relaxed structure for Si defect at the 
edge of the nanoribbon. We find that the optimization 
is accomplished mainly in two steps. Firstly, the single 
Si atom will be rejected outside the graphene nanorib- 
bon and carbon atoms around the defect undergoes re- 
construction to form a pentagon similar as edge vacancy 
defect in panel (b). Secondly, the lonely Si atom seeks 
a stable position through forming a pentagon with other 
four carbon atoms. In this sense, the edge Si defect can 
be understood as the adsorption of one Si atom to the 
relax configuration of an edge vacancy defect. So the 
edge Si defect can provide a slightly stronger scattering 
to phonon modes than the edge vacancy defect, which is 
confirmed in the transmission function shown in Fig.|3](a) 
by dot line (red online). As a result, the thermal conduc- 
tance is reduced a bit more by edge Si defect as shown in 
Fig.|4](b). Similarly as the edge vacancy defect, one edge 
C-C bond was destroyed by edge Si defect. As a result, 
the edge mode localized around this bond will show red- 
shift by the reconstruction of the edge. Fig. [1] (f) shows 
the relaxed structure for a Si defect inner the graphene 
nanoribbon. No curvature is generated by this point de- 
fect. The three bonds around the Si atom are elongated 
to be 1.7 A. The size of the Si defective region is about 
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FIG. 6: (Color online) The COOP in nanoribbons. (a). 
COOP of four electron orbitals between atom A and its three 
neighboring atoms in pure nanoribbon. (b). COOP of four 
electron orbitals between atom A in the center of vacancy 
defect and four neighboring atoms. Insets show the relaxed 
configurations. 



3.5 A. Fig.|4](b) reveals an interesting fact that the ther- 



mal conductance is not sensitive to the position of the 
Si defect. Almost a same value is obtained for thermal 
conductance of nanoribbon with edge and inner Si de- 
fect. It is because the silicon and carbon elements are in 
the same group with four valence electrons, thus no cur- 
vature is generated around Si defect. In the meantime, 
the Si defect only affects a limited range (within 3.5 A). 
Similar as the inner vacancy defect, the inner Si defect 
almost has no effect on the edge modes. 



IV. CONCLUSIONS 

To conclude, we study the phonon thermal conduc- 
tance in graphene nanoribbon by combination of ab ini- 
tio SIESTA package and NEGF approach. We suggest a 
correction to the force constant matrix to solve the con- 
flict between long-range character of first-principle calcu- 
lation and FNN property of NEGF. We first find the re- 
laxed structure of the nanoribbon and vibrational phonon 
modes. Special attention is paid to the localized edge 
modes. We then analyze vacancy and Si point defects 
in the nanoribbon and find that that the edge vacancy 
defect has limited effect on the thermal conductance as 
it only causes a further reconstruction of the edge. The 
inner vacancy defect gives the strongest scattering for all 
phonon modes by forming a saddle-like surface around 
the defect, thus greatly reduces the thermal conductance. 
For the Si defect, the thermal conductance is almost the 
same for edge and inner defect. We find that the edge 
defects can introduce a strong down-shift to the phonon 
edge mode localized around this defective bond, since the 
stronger edge bond was destroyed by edge defects. 
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